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A FINITE ELEMENT FORMULATION FOR EVALUATION OF 
CRACK BLUNTING EFFECTS IN EUSTO-PLASTIC SOLIDS 


J. R. Oslas 
INTRODUCTION 

The successful application of linear elastic fracture mechanics (LEFM) 
to prediction of fracture in metals is founded upon the representation of 
geometric flaws as infinitely sharp cracks for purposes of analysis. The 
crack model and linear elastic stress analysis provide stress intensity 
factors suitable as transfer parameters relating brittle fracture condi- 
tions in specimens and structures. Confidence in the scheme derives from 
successful correlation of failure data* which in practice obtains only 
under conditions of so-called small scale yielding wherein the flow tip 
region is subjected to a high degree of elastic constraint. 

The sharp crack model allows the use of continuum analysis, linear 
elasticity, as a basis for predicting a micromechanical process, fracture, 
by providing a characterization of loading conditions affecting a very 
small volume of material. In attempting to generalize fracture mechanics 
for application under conditions of significant yielding the appropnateness 
of the sharp crack representation must be re-evaluated. In the presence 
of such yielding not only does the relevance of linear elasticity become 
questionable, but also the involvement of larger volume of material in 
the fracture process brings the sharp crack model itself into question. 

In lieu of positing microstructural fracture criteria as a basis for 
engineering correlation of fracture data, one would like to evolve a 
fracture mechanics for non-brittle conditions which, like LEFM, employs 
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continuum analysis for definition and evaluation of transfer parameters. 
Both the J-Integral [1] and COD [2] approaches to ductile fracture pre- 
diction are consequences of this philosophy. It should be noted that, 
in principle, neither of these approaches compels the use of a sharp 
crack model . 

The model described below Is Intended to provide a tool for Investi- 
gation of the behavior of initially sharp flaws which are allowed to 
blunt as applied loading Is Increased. The approach Is motivated by the 
observation that, to the extent that actual flaws have tip radii so small 
that they may be viewed as sharp cracks, blunting must occur under load 
and will Influence internal fields over a finite region near the tip. The 
hypothesis that the blunting effect is Important 1s reinforced by recent 
experience with J-Intcgral correlations [3] In which a stretch zone 1s 
accounted for in determining a critical loading for fracture. 

The case for the utility of this form of analysis is not absolute. 
Indeed the same argument, "a sharp crack must blunt," could be applied 
to cases for which LEFM Is adequate. The distinction lies in the volume 
of material Involved in the fracture process which must be larger when 
significant yielding occurs. The fonnulation presented here will pro- 
vide a means for assessing the effects and range of Influence of crack 
blunting upon internal fields near flaw tips. Should blunting be Judged 
significant fracture transfer parameters sensitive to its effects must 
be defined and correlation of fracture data attempted before credible 
Inferences can be drawn. 
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THE BLUNTING PROBLEM 

Consider the problem of a finite length sharp crack In an elasto- 
plastlc solid. We restrict our attention to two dimensional problems. 
Under Initial load application the near tip field will be dominated by 
the characteristic Inverse square root singularity of the Williams 
eigenvalue expansion [4] for the small strain elastic problem. Were a 
closed form solution to the large strain elasto-plastic problem available, 
one would expect to observe Immediate yielding near the flaw tip, and 
i $ Immediate blunting. In lieu of such a solution our Intent Is to provide 

for numerical prediction of these effects with recognition that resolution 
will be Irrevocably limited by the spatial discretization process inherent 
In any numerical procedure. The approach taken here is to employ a 
special finite element in the flaw tip region which is designed to permit 
anticipated forms of behavior and to account for the effects of geometric 
changes upon the mechanics of the problem. 

Most previous nonlinear analyses of crack problems have incorporated 
either material or geometric (blunting) nonlinearity, not both. In 
neither case has a wholly satisfactory solution been obtained. The 
achievement of any solution at all has required approximations in problem 
definition sufficient to render the results of debatable utility. Certain 
aspects of this experience have guided development of the present TOrmula- 
tion for the combined nonlinear problem. 

1. Asymptotic near field deformation theory plasticity solu- 
tions for the small strain problem have been obtained by 
Hutchinson [5] and Rice and Rosengren [6], Solution by 
formulation of a nonlinear eigenvalue problem required 
neglecting elasticity and limiting attention to power 


law plastic materials. These results confirm thj existence 
of singular deformation fields In the presence of (deforma- 
tion theory) plasticity and suggest a characteristic form 
analogous to that of the linear elastic problem. The re- 
sults are of limited use by virtue of the power law material 
restriction and suspect by their neglect of both elastic 
effects and history dependence of plastic flow. At best 
the results are appropriate to high excitation of cracks 
far removed from exterior boundaries. 

2. The limitations inherent in the Hutchinson, Rice results will 
be removed in the context of numerical solution by the 
special finite element under development by Swedlow [7,8]. 
This approach preserves a variable order deformation gradient 
singularity for arbitrary elasto-plastic solids and employs 
an incremental theory. These minimum material behavior re- 
quirements must be preserved in an analysis of the combined 
nonlinear problem. 

3. The only available analytical solution for the blunting 
problem was provided by Knowles and Sternberg [9]. They 
consider a limited class of hyperelastic materials whose 
large strain behavior is similar to that of power law deform- 
ation theory plastic materials. Their near field, high 
excitation asymptotic solution confirms persistence of 
deformation gradient singularities in the presence of 
blunting. Similar to the Hutchinson, Rice results a rela- 
tively simple characteristic singularity is identified. 



The present combined nonlinear formulation Is founded upon the pre- 
sumption that a relatively simple singular deformation gradient field 
prevails In the vicinity of a blunting crack, a presumption consistent 
with the foregoing experience. Furthermore, In view of the fact that 
a formal near field expansion solution Is not within reach, there Is 
little motivation for restrictive simplification of the material model 
beyond that undertaken to permit economic computation. The analysis 
is based upon a generalization of conventional Jg flow theory 

which Is appropriate to the large strain elasto-plastic problem for an 
arbitrary Isotropic hardening material. 

A final Important feature of the present formulation Is suggested 
by the experience of this author [10,11] and others [12] In attempting 
solution of the combined nonlinear problem without special attention 
to modeling the near field. At Issue is resolution of the effects of 
blunting upon Internal fields. In small strain finite element crack 
analysis, a special tip element may be sized in terms of crack length. 

Once an element size Is determined appropriate to a particular formula- 
tion that size may be used in a variety of different problems. In 
crack blunting analysis another size scale becomes Important, the "size" 
of the blunted crack tip. If tip size is thought of 1n terms of a tip 

Inscribed radius, as suggested by Srawley [13], the size scale of interest 

-5 -1 

may range vn'th load level from 10 to 10 times the crack length in the 

course of a single analysis. The change of size scale with load will depend 
primarily upon stress-strain curve and to a lesser extent upon loading 
type: tension, bending, etc. It Is clear that no one special element 
size can give results of equivalent accuracy over any substantial range 
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of excitation in a single analysis, much less from one problem to another. 

The problem may be treated in either of two ways. The range of Influence 
of the blunting field singularity may be taken as a variable determined 
by the analysis, or a single value may be used which is appropriate to 
the load level at which results are desired. The utility of the latter, 
"simpler," approach cannot be evaluated a priori as an unknown degree 
of history insensitivity is presumed. No unique solution to this dilenma 
has presented Itself. The formulation outlined below allows the spatial 
range of Influence of the blunting field sinqularity to be either an 
unknown, or a fixed value, without influencing the compatability charac- 
teristics of the finite element model. Evaluation! of the influence of 
this quantity upon analysis results must be of high priority once an 
operational computer program is available, 

SPECIAL ELEMENT FORMULATION 

The elasto-plastic crack blunting finite element formulated below 
is designed for implementation in a host finite element program capable 
of solving small strain elasto-plastic problems employing a J 2 flow theory. 

In principle the host program may employ elements of any order. The 
present discussion presumes use of an 8 or 9 node isoparametric quadri- 
lateral such as that being developed by Marino [14]. Coupling between 
host program and special element is in tenns of nodal variables (dis- 
placement, force) only. The special element and its associated displace- 
ment interpolation functions are defined in Lagrangian terms, i.e., on 
the undeformed geometry. Finite deformation effects are considered 
within the special element, and may be neglected elsewhere, since the special 
element employs a large strain elasto-plastic representation which reduces 
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smoothly to familiar plasticity theory in the limit of infinitesimal 
deformation. 

The following discussion is limited to development of the essential 
features of the special element. Straightforward operational detail 
prerequisite to writing a computer program is not included. 

Definition and implementation of the special element requires establish- 
ment of three distinct, consistent formulations. 

1. A field problem with an associated variational 
principle. 

2. An elasto-plastic constitutive model for large 
strains. 

3. An interpolation model for the special element. 

Definition of a special element interpolation model incorporating key features 
of the near crack tip region necessitates use of a Lagrangian coordinate 
frame. It also requires introduction of algebraic quantities whose evalua- 
tion necessitates availability of a variational principle. Consequently, 
the large strain formulation developed by Osias [15] and Osias and Swedlow 

[16] is inappropriate to items 1 and 2 above. An alternative, and essen- 
tially equivalent, Lagrangian formulation has been developed by Hutchinson 

[17] and will be employed here. The formulations are discussed and com- 
pared by Key, Osias, Belytschko and Hutchinson [18]. 

Field problem : We choose a convected coordinate model as defined by 

Green and Zerna [19] and Nemat-Nasser [20], The' initial geometry, interior 

★ I 

Bjj and boundary 9B^, is described by material coordinates X in an ortho- 
gonal frame with metric 6jj and base vector Gj. The deformed configura- 
tion, B bounded by 3B, is described by coordinates in a convected frame 

^Conventional indicial notation is employed with all indices varying over 
1,2,3. The summation convention is employed. Partial differentiation is 
indicated by a comma and covariant differentiation by a semicolon. 
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with n»tr1c and base vectors g^. For a material particle the dis- 
placement field 1s 

IT • x^g^ - X^Gj - U^Gj 

In (1) we choose to Identify the material frame tensor components as will 
be the case throughout this development. Total strain 1s characterized 
by Green's strain (2). 

ElJ • r <«I,J " 

Of the several possibilities It Is convenient for our purposes to 
employ the so-called KIrchoff stresses 


Cl) 


( 2 ) 


•IJ - fzm _I0 


e i/gTS 


(3) 


IJ 


(4) 


In (3), g = |g^j| and G = ^ material frame components 

of Cauchy stress. Thus for a surface whose undeformed state outer normal 
Is N = N Gj the deformed state traction is found as: 

T = T^Gj = 

Restricting attention to the case of static problems without body forces 
equilibrium equations may be established in the form (5). 

(s'^’ + s“-u!l).|^ = 0 (5) 

Large strain elasto-plastic quasi-static behavior is governed by 
rate equations which are first order homogeneous in time. These equations, 
of the form (6), motivate formulation of an analysis as an initial- and 
boundary-value problem for rate quantities. This problem, also first 
order homogeneous in time, is quasi-linear and, as for the small strain 
problem [21,22], provides a basis for efficient incremental solution. 


•IJ ^ pIJKL^ 


KL 


( 6 ) 
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In (6) S are the components of the convected rate of Kirohoff stress 
and E|j are components of the material derivative of Green's strain. As 
will be shown shortly, the constitutive tensor P In (6) possesses the 
symmetry: 

pIJKL . pKLIJ (7) 

The rate problem Is formed employing rate equilibrium equations 
(8) and boundary conditions (9) for the quasi -static problem. 

( 8 ) 

Equilibrium equations (8) apply to the undeformed Interior configuration 
and boundary conditions (9) to the undeformed boundary 38^. are 
the velocity-field components. 

A basis for finite element solution may be obt'.lned by establishing 
a variational statement of the problem defined by equations (6-9). Straight- 
forward manipulation yields the stationary principle (10) for the case where 

•T 

T Is prescribed on 


«I « 0 


I 


I 


; . vkui „ ,h„ 



( 10 ) 


In (10) the variation Is taken with respect to the velocity field, noting 
that Green's strain rate Is given by: 


IJ 


1 

J 


(V 


I;J * '^J;I 


V 

^;rK;J 


^;0''k;I^ 


( 11 ) 
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Finite element analysis employing this variational formulation requires 
definition of a large strain consitutive model providing the symmetry (7) 
and an Interpolation model for the material frame components of velocity 
within the elen^nt. The analysis will perforce take the form of an In- 
cremental procedure Involving evaluation of the velocity field at a 
sequence of times during the course of the deformation process. 

Elasto-plastic Constitutive Equations ; The constitutive model 
which we will employ preserves all of the features of conventional Iso- 
tropic Og theory of plasticity. The model utilizes an arbitrary 

work-hardening stress-strain curve developed from tensile test data and 
Includes provision for treating elastic unloading. Similar to the work 
of Hibbitt et al. [23] and Os las [15]. objectivity 1s preserved by Intro- 
duction of a Jaumann stress rate. The formulation Is restricted to small 
elastic deformations. 

The model is developed by Hutchinson in [17] and the derivation 
will not be repeated In detail here. Elasto-plastic flow equatloh* are 
established by decomposition of the Green's strain rate Into elastic and 
plastic components. Elastic deformation Is described by an isotropic 
hypoelastic generalization of Hookean elasticity. Yielding and plastic 
flow are controlled by a loading function taken as the second deviatorlc 
Invariant of the Kirchoff stress (3). 

A generalized expression for total Green's strain rate Is found In 
the form: 

• 1 *KL 

Mj ■ t 


( 12 ) 
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where a ■ 1 for plastic flow and a * 0 for elastic behavior. In (12) 
(E,v) are the constants of linear elasticity, gjj^ Is the metric of the 
convected system x^,s*^ are the devlatoric components of and f Is 
a work hardening parameter taken as a function solely of Jg* 

->2 * (J-'SikSjL*"*”’ 


(13) 


1 . n « 

^2 9ik9jl* ^ 


KL 


In (12) and (13) the Oaumann stress rate S p»*eserves flow rule objec- 


tivity and Is related to the convected rate by (14) 


jlj . jlO ^ glV4. 


■f 


'KL « ^ 'KL 

Substitution of (14) In (12) and subsequent Inversion obtains constitutive 
equations of the form (6) and provides the symmetry requisite for the 
variational form (10) to hold. 


(14) 


IJ_KL 


pIJKL . pKLIO . ^ (glyL , gILgJK, . ^ 


(13/ 


a _IJ,KL-i 
-gs s J 


- (^)[9'VL + + g"^S« * 

Material property Information enters (12) and (15) through the elastic 

— ^ 

constants, which take values ascribing to the undeformed state, and through 
the work hardening functions f and G. Evaluation of (12,15) for simple 
tension yields the relationships (16) 

[1 + j • «5 - ‘ 

- (|)’'^(r * I u - 2^1 


( 16 ) 
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In (16) is the instantaneous slope of the curve relating true stress o 
and logarithmic tensile strain e. v = is the instaneous contrac- 

tion ratio, e 2 being logarithmic transverse strain. 

Special Element Interpolation Model : The model is developed tc pro- 
vide a number of kinemati features judged to be necessary for representa- 
tion of crack blunting under Mode I loading. It is intended for use with 
the variational formulation (10) and does not introduce kinematic incom- 
patibilities anywhere in the domain of the analysis. 

The near crack tip total displacement field is represented as a 
Combination of functions incorporating nonsingular, and variable order 
singular, gradients. Material frame components of displacement are 
approximated as functions of the undeformed coordinates so as to pre- 
serve reference to the original crack geometry. The singular gradient 
portion of the displacement iwdel is further characterized by introduction of 
an explicit length scale which may either be held at a fixed value or 
allowed to vary during the analysis. This parameter is intended to allow 
the range of influence of the singular gradient field to vary as crack 
blunting occurs. 

A velocity field approximation appropriate to formulation of a rate, 
or incremental problem is obtained by differentiation of the displacement 
model with respect to time. Whereas some of the algebraic variables- in 
the displacement model appear in a nonlinear fashion, their time rates 
appear linearly as unknown coefficients in the velocity representation. 

Thus a basis is provided for formulation of an incremental finite element 
analysis requiring solution of linear equations. 

The special element domain is defined on the undeformed state and 
is the region 0 ^ e ^ ir, 0 ^ R/R^ < 1 of Fig. 1. (R,e) are undefomed 


material, or Lagranglan, coordinates and Is fixed. The value of 

must be large enough to Insure that the region within which singular 

behavior Is affected by blunting Is contained within the special element. 

.9 

In application R will be of order 10 a, where a Is crack length. The 
domain Is partitioned Into N pie-shaped subregions, or base elements. 
Figure 1 shows the case N = 3; actual application will require N to be 
of order 10. 

The total displacement field Is approximated as: 

u,(x’.x^.t) - Ujj(x’.x^.t) + Ujj(R.e.t) (18) 

* 

where U^| are defined independently within each base element and have 
nonsingular gradients. The are defined over the full special element 
domain and contain terms leading to singular gradients. 

For null U^j the special element would consist of an array of con- 
ventional higher order, nonsingular finite elements surrounding the 
crack tip. The displacement Interpolation functions U^j within each 
such element are constrained only by the requirements of geometric and 
displacement compatibility with each other and with host program elements 
across R/R^ = 1.0. Figure 1 suggests the use of 6 node isoparametric 
triangles; the base elements could also be degenerate, nonsingular, 8 or 
9 node isoparametric quadrilaterals. The “oI must account for rigid 

motions and uniform deformation states. Without prejudice as to specific 

th * 

form we will denote these functions for the m base element as 

l/”,(x\x^,t) = i A^^(t)(|>j (x\x^) (no sum on I) (ig) 
a=l 

★ — — — 

Henceforth Latin indices o,s,m are reserved for use as defined here and 
do not denote tensor components. 
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where time dependence appears through unknown coefficients A"jt). and 
n Is the number of nodes associated with the base elen«nt. 

The special el&nent displacement interpolation model is completed 
by specificatim of “si which is defined over the domain 0 £ R/R^ £ 1.0, 

0 ^ e ^ IT, where Rj Is to be determined as part of the analysis. 

For R ^ R^, and within base element m, is written as: 

Usi(R,e) = (20) 

for Rj/Rell-O* R/Rgll-O 

(li - 1 )ae <_ 0 mA0, 0 <. ^ a 0 ae = n/N 

The precise functional form of the functions Rj and 0j Is somewhat 
arbitrary. If an exact asymptotic solution for the blunting problem were 
available it would be employed here. In lieu of such a solution we 
employ (20) with the functional forms (21-23). 


R, = 


n.(0) 

(|-) [1 




( 21 ) 


where 


nj(e) = + n^j sin 0 + ngj sin 2e + ... + Opj sin pe (22) 

I I 

0?^Om) = ^ Bj sin Y0„ + Z C? cos y0„ (23) 

I m 1 Iy 1 Iy hi 

Y=| ' y=l . ' 

where the b'J^ and Cj^ are time dependent. 

Thus complete definition of the model requires in addition to specification 
of N and Rg, the choice of lengths for the series (22) and (23). The 
nunter of unknown time dependent coefficients in (23) is reduced by re- 


quiring that (24,25) hold at the radial boundaries between base elements. 

e" = ef ’ (24) 




15 


36 

m 




36 , 


m+1 


(25) 


The form of the displacement field model U^j of (20) as detailed In 
(21-25) f**llows from provision for a variable order displacement gradient 
singularity in the vicinity of the crack tip. The singularity Is Introduced 
such that the special element displacement field is single-valued and 
further such that at R/Rj = 1.0 for 0 ~ 


’“si . 

TR 


3e -“sI-O 


(26) 


1,e., the singular gradient field vanishes smoothly as R which may 
vary during the analysis. 

The velocity field interpolation model for use in (10) is established 
by differentiating (18) and provides for each component at time t: 

Vj(x',t) = U^J + Ujj (27) 

As a consequence of the form (20) of the displacement model, the velocity 
field will be single valued and will couple smoothly to the host program 
velocity model across R/Rg = 1.0. Furthermore, the singular gradients 
terms, U^j and its first derivatives, vanish at R/Rg = 1.0. The velocity 
field interpolation model may be written within the m ^ base element as 

(no sum on I) (28) 

(29) 


V, . ujj . Rjef . 


3R 


I 


3 Rt p 3n.(e) . 3 Rt . 

^ (no sum on I) 

I an ol an an^ 6l aRg s '' ' 


©'J = z Bf sin Y0 + z C? cos ye 
I Iy m 1 lY m 

Y-1 Y=1 


0=1 


(30) 


( 31 ) 


Thus at any time during the deformation process the deformed geometry 
is established through the displacement field (20) by coefficients Aj^ 
of (19), Hgj of (22) and Cj^ of (23) as well as the value of 
appearing in (21). Presuming these quantities to have been established 
by preceeding analysis, the velocity field is found by substitution of (27) 
Into the functional (10) which will be quadratic in the time rate quanti- 
ties appearing as coefficients in (28-31). The variational processs will 

*m *Ri * * 

provide linear equations for (Aj^» Bj^, Cj^, R^, n^j, ngj), the unknown 
coefficients. 

Solution for the time rates and integration over an increment of 
time will establish new values for the displacement field parameters, 
allowing a new rate problem to be defined at the later time. 

In this manner a complete deformation history may be established by solu- 
tion of a sequence of linear algebraic problems. 

BLUNTING MODEL APPLICATION 

The above formulation will permit solution of crack blunting problems 
within the framework of conventional incremental finite element analysis. 
Computation cost will be directly dependent upon the number of base 
elements anployed within the special element as well as the number of 
terms carried in the series portions of the interpolation model. The 
incremental stiffness equations will be linear and well -conditioned. 
Solution may employ conventional reduction techniques. 

Problem definition will require specification of loading history, 
geometry, and material properties as well as non-zero initial values of 
the coefficients in nj(e) (22) and Rg in (21). Typically one would employ 
the small strain elastic values for the coefficients, i.e., n^j = 1/2, 


Hgj = 0 (0 * 1,..., p). An Initial value for Is less obvious. A 
value of the order of the osculating circle for the deformed elastic tip 
radius seems appropriate. The particular value giving the most accurate 
elastic stress Intensity factor would be a viable starting point. Solu- 
tion sensitivity to initial values of will have to be investigated by 
numerical experimentation. Should variability of prove to be an unr 
desirable feature the model permits Its value to be fixed. 

The special element will permit blunting problem solution within the 
constraints Implied by the form of the Interpolation model. The ad hoc 
form of this model makes it most appropriate for evaluation of integral 
parameters, e.g., crack tip displacements or J-integral values. Results 
for internal field quantities will be directly controlled by the form of 
the particular interpolation functions employed. The utility of the analysis 
must ultimately be evaluated on the basis of its ability to support correla- 
tion of fracture data for ductile materials. 
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FIGURE i: SUBPARAMETRIC ELEMENT CONFIGURATION 
FOR CRACK BLUNTING ANALYSIS 
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